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New quantum-mechanical phenomenon in a model of electron-electron interaction in 

graphene. 
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Budker Institute of Nuclear Physics of SB RAS, Novosibirsk, 630090 Russia and 
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A quantum mechanical model of two interacting electrons in graphene is considered. We concen- 
trate on the case of zero total momentum of the pair. We show that the dynamics of the system 
is very unusual. Both stationary and time-dependent problems are considered. It is shown that 
the complete set of the wave functions with definite energy includes the new functions, previously 
overlooked. The time evolution of the wave packet, corresponding to the scattering problem setup, 
leads to the appearance of the localized state at large time. The asymptotics of this state is found 
analytically. We obtain the lower bound of the life time of this state, which is connected with the 
breakdown of the continuous model on the lattice scale. The estimate of this bound gives one a 
hope to observe the localized states in the experiment. 

PACS numbers: 73.20.Mf, 73.22.Pr, 03.65.Gc, 03.65.Nk 



I. INTRODUCTION 

Nowadays a great deal of effort has been devoted to the 
experimental investigation of the transport properties of 
graphene, see recent revievi?^. One of the important re- 
sults of these experiments is the observation of high mo- 
bility of the charge carriers^. Many papers have been de- 
voted to the investigation of the influence of the electron- 
impurity interaction on the mobility of the charge carri- 
ers, see reviews^i^. Study of the electron-electron inter- 
action in graphene is also important for understanding of 
this effect, see ReL-. However, a theoretical progress in 
this problem is rather limited^.'.!. 

It is well established now that the low-energy single 
electron dynamics in graphene is described by a massless 
two-component Dirac equation^ii^""— 

ihdtip {t, r) = hip (i, r) , 

where the hamiltonian h has the form 

h = vpcr ■ p, 

vp is the Fermi velocity, p — —ih^, and cr — {ox,(Jy) 
are the Pauli matrices acting on the pseudospin variables. 
Below we set h = vp = 1- Evidently, the pair of non- 
interacting electrons can be described by the equation 

idtip {ri,V2,t) = Hoip{ri,r2,t) , (1) 

Ho^hi + h2 = o-i • pi + 0-2 • P2 , (2) 

where ip (ji ,r2,t) is the wave function of the system, de- 
pending on the coordinates and pseudospin variables of 
both electrons. The generalization of Eq. ([Ij to the 
case of interacting electrons is a highly nontrivial prob- 
lem. The origin of the difficulties is the necessity to 
take into account the interaction with the electrons be- 
low Fermi surface. This interaction results in the ex- 
istence of the electron-hole excitations in the interme- 
diate states. The account of the corresponding effects 



in quantum electrodynamics (QED) leads to the Dyson- 
Schwinger equation (which, for the bound states, reduces 
to the Bethe-Salpeter equation), see, e.g., RefJ^. How- 
ever, in the nonrelativistic QED systems, the effect of 
virtual electron-positron pair in the intermediate states 
is small. For massless electrons in graphene, the nonrela- 
tivistic approximation is not applicable and the effect of 
virtual electron-hole excitation may be crucially impor- 
tant for the problem of electron-electron interaction. The 
approach based on the Bethe-Salpeter equation was used 
in Ref jii in the investigation of electron-hole interaction 
in graphene. 

Though the influence of the electron-hole excitations 
can be very important, nevertheless, as a first step 
in the investigation of the electron-electron interaction, 
it makes sense to ignore this effect and to model the 
electron-electron interaction by replacing Hq — > Hy in 
Eq. dD), where 



m 



Ho + V{r)=(Ti-pi+(T2-P2 + V (r) (3) 



and V (r) — V{\ri — r2\) is the electron-electron inter- 
action potential. Recently, this model has been consid- 
ered in RefJ^, where the eigenfunctions of Hy have been 
analysed. The solutions found in Ref.^^ appeared to have 
unusual properties. In order to understand the origin of 
these properties, we revisit in the present paper the so- 
lution of the stationary equation Hyip = Etp. We also 
consider the time-dependent problem and demonstrate 
that the unusual properties of the eigenfunctions of the 
hamiltonian are reflected in the very specific properties 
of the time evolution of the wave packets. 

We restrict our consideration to the specific case of zero 
total momentum of the pair and search for the solutions 
being the eigenstates of the operator 



J^ 



{al+(jl)~id^ 



where Lp is the azimuth angle of the vector r = ri — r2. 
We assume that the potential V (r) is a smooth positive 



monotonically decreasing vanishing function. To include 
the important case of the Coulomb potential into the con- 
sideration, we allow for the r~^ growth of the potential 
at r — >■ 0. The solution of the stationary equation shows 
that the wave functions with J^ = are smooth func- 
tions for any energy E. This is also valid for the wave 
functions with Jj, ^ and the energy above the maxi- 
mum of the potential Knax = V (0) or below zero. For 
Jz "/= and < -E < Knax, the solution of the stationary 
equation necessarily has singularity at the point r-^ (E) 
determined by the condition 



E = V{r,), 



(4) 



which is in agreement with Ref.^^. Such a behaviour con- 
tradicts a common wisdom which tells one that the wave 
function should be a smooth function in the region where 
the potential is also smooth. We show that the existence 
of the singularity in the wave function is related to the 
degeneracy of the derivative matrix in the hamiltonian. 
We find an important new feature of the energy spec- 
trum: the additional degeneracy of the states with fixed 
J^ ^ and energy in the interval (0, Vmax)- 

For the time-dependent problem we choose the ini- 
tial conditions corresponding to the wide spherical wave 
packet with fixed J^ and the average energy Eq (average 
value of the hamiltonian) moving toward the origin from 
the large distance rg ^ A, where A is the width of the 
packet (the energy dispersion in the packet ~ A"-'^). The 
direct numerical calculation reveals a remarkable picture. 
At rather large time t > ro one observes not only a re- 
flected wave packet moving toward the large r, but also 
a narrow peak in the vicinity of r^, {Eq) with the width 
oc 1/A. The total norm of the wave function is conserved 
as it should be for a hermitian hamiltonian. In order to 
check consistency of the results obtained, we demonstrate 
that the solution of the time-dependent equation based 
on the decomposition of the initial wave packet over the 
stationary wave functions reproduces the direct numeri- 
cal solution of this equation. Using this decomposition, 
we find the large-time asymptotics of the emerged peak 
analytically. 



II. GENERAL PROPERTIES OF THE MODEL 

Obviously, the hamiltonian H in Eq. (|3]) commutes 
with the total momentum P = Pi + P2, and we can 
search the wave function in the form ?/'(ri,r2,i) = 
exp (iPp • R) ?A (i, r) , where R = (ri-f-r2)/2 is the 
center-of-energy coordinate and r = ri — r2 is the rela- 
tive position vector. Note that the wave function ip {t, r) 
depends nontrivially on the system total momentum Pq, 
see Refii^. Below we consider a specific case Pq = 0. 
Then the wave equation has the form 



H^{rTi-rT2)-p + V (r) 



(5) 
(6) 



where p = — iV. The hamiltonian H commutes with the 
operator 

j'^S^ + L^ = ^{at + a^^)-td^ (7) 

(if is the azimuthal angle of the vector r) and the operator 



O = S^ - 2 



{S')\ 



(8) 



where S^ = i Z^Li {'^i + '^2) ■ The operators J^ and O 
also commute with each other. Therefore, we can search 
for the solution of Eq. ([Sj to be the eigenfunction of J^ 
and O: 

^o{t,r) = e^*^^(aoo(i,r) |0, 0) + e-^'^an (i,r) |1,1) 

+ e«^ai_i(i,r)|l,-l)), (9) 

V'2(t,r) = e'*%(t,r)|l,0) (10) 

so that J'^ipk = Mtpk and Ot/jk ~ kipk- Here |s, Sz) is the 
eigenfunction of the operators S^ and 5*^. It is convenient 
to pass from the functions aij to the functions 

.aii+ai_i .flu— ai_i 
/=* -1^ , h = i -= , d = aoo. (11) 

Using Eqs. ©, ©, dTU]), and dm, we obtain 

idtg^V{r)g, (12) 

(13) 
idth = V{r)h~ 2drd, (14) 

idtd = V{r)d~—f + 2(dr + -]h, (15) 
r \ r J 

The last three equations can be represented in the matrix 
form 



idtf = V{r)f~—d, 
r 



idtF = HrF, 



(16) 



where 



F= [h] , Hr^ { V{r) -29^ . (17) 

\d) V-^ 2(9. + i) Vir)J 

It is easy to see that Hr is a hermitian operator, i.e. 

00 00 



drrFlHrF2 = / drr 



Fo 



for continuous functions Fi^2 [r), decreasing sufficiently 
fast when r — >■ 00 and finite at r — Q. 
The general solution of Eq. ([12]) is 

(18) 



ff(t,r)=.g(0,r)e-^^M*, 



whereas the general solution of Eq. (IT51) can not be found 
analytically. 



Conserved current and density. The conserved cur- 
rent and density for Eq. ([5]) have the form 

j=^jH(7l-CT2)ij, p = i>H. (19) 

For two sohitions ip^ and ^1^2, Eq. (jlOl) . the current and 
density are expressed as 

^o: > = 4Im(d/i*), i^ = -4Re(dr), 

p=\f\^ + \h\^ + \d\\ (20) 

V'2: > = 0, j^ = 0, p-|5l', (21) 

where jr and ji^ are the radial and angular components 
of the current, respectively. 



III. STATIONARY PROBLEM 

Let us first consider the stationary equation for the 
function g: 



Eg^V(r)g. 



(22) 



This simple consideration helps one to understand better 
the properties of the solutions of the stationary equation 
for the function F. For E > VJnax, the equation (1^^ has 
no solutions, while for < -E < VJnax its formal solution is 
9a (r) = ^ (r — a), where a is determined by the equation 
E — V (a). The functions ga (r) for different values of a 
are mutually orthogonal and normalized by the condition 



dr r ga (r) ga (r) = aS{a- a) 



(23) 



Note that the density p{r) = \ga (r)| , Eq. (pTj) . corre- 
sponding to this solution, is not well-defined. Neverthe- 
less, the functions ga (r) form a complete set and can be 
used to solve the time-dependent problem. Indeed, 

oc 

g (i, r) = [da e-'^'^^^'g (0, a) ga (r) = g (0, r) e^^^^'')* 



in agreement with Eq. (|18p . 

Let us now pass to the consideration of the stationary 
equation 

(24) 



EF = HrF . 



The hamiltonian Hr is a first-order differential operator, 
see Eq. (1171) . It is known from the theory of ordinary 
differential equations that the solution y of the system 
dry (r) — A(r) y (r) can have singularities only in the 
points where the elements of the matrix A (r) are singu- 
lar. We can not, however, represent Eq. (|24p in this form 

/O \ 
since the matrix —2 in front of the derivative dr 

\0 2 J 

in Hr is degenerate. We show below that this degeneracy 
leads, for M 7^ 0, to the appearance of the singularity of 
the solution F in the point r = r^,. 



Second- order equation for d. The explicit form of Eq. 
([M]) reads 



iE^V)f 



2M 



{E-V)h^ -2drd, 

2 2Af 

(E-V)d^ 2drh +-h / . 

r r 



(25) 
(26) 
(27) 



Using the first two equations in order to eliminate the 
functions / and h from the last equation, we obtain 



d" +p{r)d' + q{r)d = 0, 
V 1 

E -V^ r' 



p{r) = 



q{r)^-{E-Vf-^' 
4 r^ 



(28) 
(29) 

(30) 



where a prime denotes the derivative with respect to r. 

Boundary condition at r = 0. Let us determine the 
boundary condition at r = 0. If V (0) < 00, the general 
solution of Eq. (|^5|) behaves near r = as 



—'•"--;::'• :;::^ ™ 



while for the case of Coulomb singularity, when V (r) — > 
a/r, the asymptotics of the general solution has the form 

dKair''-^/^ + a2r-''-^/^ , (32) 



where u = ^y/AM"^ -f- 1 — a^. Here ai and 02 are some 
constants. We choose the boundary condition at r = 
as 



02=0. 



(33) 



This condition provides that ^„drr\F\ < 00 for suffi- 
ciently small (5 > 0. 

Analytical properties. To understand the properties 
of the solution d of Eq. (|28|) , we consider the analytical 
properties of the coefficients p{r) and q(r), Eqs. (|29l) . 
([501) . on the interval [0,cx)). In the origin, the coeffi- 
cients behave as p{r) ^ r"^, q{r) ~ ?'~^, so that the 
point r = is a regular singular point of the differential 
equation ([^5]) . The coefficient q (r) tends to a constant 
when r — >■ CXI, therefore the point r = cx) is an irregular 
singular point of the equation. The above properties of 
Eq. (P5)l (singularities of p (r) and q (r) at r = 0, 00 and 
boundary condition for d (r) at r = 0) are analogous to 
those of radial Schrodinger equation. The new property 
of Eq. (l28l) is the singularity of the coefficient p (r) at 
r = r^,[E), see Eq. (g]), when E £ (0,Knax)- In prin- 
ciple, the singularities of the coefficients of the equation 
do not necessarily lead to the singularity of the solution 
(and its derivatives). One can check that the general so- 
lution for M = is, indeed, a smooth function at r = r^,. 
Therefore, we concentrate on the case M 7^ 0. The gen- 
eral solution of Eq. (pS)) in this case is not smooth at 



r = r^,, which can be readily seen from the asymptotics 
of the solution in the vicinity of r^: 



d{r)^h{l + 



AP (1 - r/r,f 



ln|l - r/r*| 



+ 62(l-r/rO' 



(34) 



Therefore, we search for the solution separately in two 
regions 

d (r) = J^i^'" '^^^ "*" ^2C?rog (r) , < r < r^ 

where bi^2 and &1.2 are some constants, and the functions 
din- and dreg have the asymptotics 



dirr (r) 



4og (r) 



M2 (r* - rf 



2rl 



(r* - r) 



ln|l-r/r^|, (36) 
(37) 



at r -> r^,. 

Matching conditions at r = r^,. In order to determine 
the general form of the solution of Eq. (P5|) . we need 
to apply matching conditions at r = r^,. Conventional 
analysis of Eq. ([25]) in the vicinity of r* leads to the 
requirement of the continuity of d and d' aX r — r^, . Using 



J 



these conditions we end up with 
&i = 61 . 



(38) 



Note that 62 remains a free parameter, which means that 
the function d (r) in the region r > r* is not entirely 
determined by that in the region r < r^,. 

The boundary condition at the origin, Eq. ((33)) . fixes 
the ratio 



62/^1= /3, 



(39) 



where /3 is the constant which depends on the energy 
and the form of the potential. Therefore, we have two 
conditions, (pS]) and (|39)) . for four constants 61,2 and 61, 2- 
It means that for any energy in the interval E € (0, Knax) 
there are two linearly independent solutions which we 
choose as 

di{r) = 5i Kt (r) + /?drcg (r)] , (40) 

d2{r) = 62^ (r - rO dreg (r) , (41) 

where 9 (x) is the Heaviside step function. 

Solutions of the system (|25p - (|27p . Let us now return 
to the initial system, Eqs. ([25]) -([271). Note that Eqs. 
(P5)) and (|26p determine functions / and h up to the 
generalized function localized at r = r^,. Substituting 
(|40l) and (gl]) in Eqs. ([251) -((271), we obtain two solutions 
of the equation ([24]) : 



2M 



FiiE,r) 



F2{E,r) 



dijE.r) 

r ^ V(r) 

2drdi{E,r) 

V(r)-E 

di {E, r) 

2Md2{Ejr) _ 2rd^d2(E,r^+0) 
r{V{r)-E} M 

2d^d2(E,r) 
V{r)-E 

d2(E,r) 



5{V{r)-E) 



(42) 



(43) 



r 



Here P - stands for the principal value defined as 



1 



V{r)-E 



1 



1 



V {r)-E + iQ V{r)- E-iO 



In order to check that Fi and F2 are the solutions of Eq. 
()24l) in the vicinity of r* , one can integrate the equations 
(P5)) - (P7|) over r from r^ — 61 to r^, + (S2, and consider the 
limit (5i,2 -J' +0. 

Similar to the solutions ga (r) of the equation ([22]) , the 
functions Fi^2 contain generalized functions. The density 
p ~ F^F, corresponding to the solution Fi is not inte- 
grable at the point r = r^,, while that, corresponding to 
F2, is ill-defined. Therefore, for energies in the interval 
(0, Knax) there is no solution F {E, r) of the stationary 
equation with the density being an integrable function in 



the vicinity of i\ (E) . This statement is in clear contra- 
diction with the statement of Ref J^, where it was claimed 
that nonanalyticites at r = r* give a finite contribution 
to the probability. 

It may seem that the same consideration of the case 
M — will also lead to the two-fold degeneracy of the 
spectrum for < E < Knax- However, it turns out that 
the substitution of ^2 from Eq. (|^T|) to the original sys- 
tem, Eqs. (p5)) - (|27)) . leads, for M — 0, to appearance of 
the ^-function term violating Eq. (I?7l) . 

For completeness, let us also discuss the properties of 
the solution of Eq. ((24|) for the energies above the max- 
imum of the potential or below zero. In this case there 
is no singularity in the coefficient p (r) on the interval 
(0, 00). The energy spectrum is not degenerate for E < 



and E > Knax since the solution is defined uniquely (up 
to the normalization) by the boundary condition at the 
origin. This solution has the form (|42p where one can 
omit the P symbol. In what follows we assume that 
Fi (E, r) ioT E < and E > Vniax is normalized as 



/>oo 

/ drrFl{E,r)Fi{E',r)^2nS{E~ E') , 
Jo 



(44) 



so that the large-r asymptotics of Fi {E, r) has the form 



Fi{E,r)''^^ 




(45) 



at e — >■ ±0. The limit depends on the sign of e (see below) 
and we denote the corresponding solutions by the lower 
index + or — , respectively. The equation for the function 
d (r) has the form (pS)) with the replacement i? — >■ E + ie. 
For e 7^ 0, the coefhcient p (r) of this equation has no 
singularities on the interval (0, oo) and its solution d± is 
fixed, up to a constant factor, by the boundary condition 
at the origin. The first two components of the functions 
F± can be expressed in terms of d± so that 



F±{E,r) 



f±iE,r) 
h±{E,r) 
d±{E,r) 



2Md±{E,r) 

r(V{r)-E^iO) 

2d^d±(E,r) 

V{r)-E^iO 

d±(E,r) 



(47) 



where (f is some function of the energy. 

Alternative derivation. We present now an alterna- 
tive derivation of the solutions (gl]) and (gS]) for M ^ 0, 
which allows one to understand better the appearance of 
the second solution F2. For this purpose we interpret Eq. 
for real E as a limit of the equation 



[E- 



HrF 



(46) 



where we assume that the limit e — > ±0 is already per- 
formed. On the interval {0,r^,{E)) the functions F+ 
and F- coincide. Without loss of generality, we can 
choose them to be real on this interval. Then, obviously, 
F+ {E,r) = FZ {E,r) on the whole interval (0, 00). On 
the interval {r^, (E) , 00) the functions F± gain imaginary 
parts which is clearly seen from the asymptotics of the 
functions d± in the vicinity of r^,: 



The prescription ±iO determines the choice of the logarithm branch for r > r*, so that 



d± {E, r) « &i 1 



AP (1 - r/r,f 



In 



r*-r 



ri. 



PTi- 



[r - r*) (1 - r/ri,Y 



(48) 



(49) 



r 



Using this formula, as well as the identity 



1 



= P 



1 



±iTT5{V{r)~E) , (50) 



V{r)~ETiO V{r)-E 

one can check that the real part of F± is proportional to 
the function Fi, Eq. (g^ . and the imaginary part of F± 
is proportional to the function F2, Eq. P5)) . 

Orthonormality and the dual basis. For the potential 
decreasing faster than 1/r, the asymptotics of the func- 
tions d± at r — >■ cx) has the form 



d± (E, r) 



2V^ 



(^^T^{Er+^)/2^^^±r(Er+^)/2^ ^ (5^) 



where c, 7, and ip are, in general, some real-valued func- 
tions of the energy. For the potential decreasing at 
r -^> 00 as a/r, one should perform the replacement 
tp ^ if — ahvEr. We choose the overall normalization 
constant c to be equal to unity, so that the asymptotics 
of the functions F± has the form 



F±{E,r) 



2^ 



±1 I eTHB'-+^)/2 + ^ 




^±i(Er+ip)/2 



(52) 



Then the functions F± satisfy the relation 

oo 

f dr rF^ {E, r) F„. (E' , r) = 2tt6 {E - E') N,,, , (53) 



N++ = iV__ 



1+72 47rM2|6i|^ 

N+- = iV_+ = 7 , 



(54) 
(55) 



where 61 = d±(E,ri,), see Eq. (j49|) . There is one 
subtle relation between constants 61 and 7, which fol- 
lows from the conservation of the total radial current 
Jr {r) — 2TTrjr, where jV is defined in Eq. pO|. Namely, 
using the equality Jr (r — ?> Tj, + 0) = Jr (r — >■ cxd), we ob- 
tain 



so that 



AnM^lhl'^ _ 1-7^ 
r*\V'{r,)\ ^ 2 



Ar++ = Af__ = 1 . 



(56) 



(57) 



Note that J^ 7^ in the region r > Vi,. The existence of 
the solutions with the nonzero radial current is the con- 
sequence of the spectrum degeneracy. For r < r^we have 
Jr- = 0. It may seem that such a behaviour contradicts 
the continuity equation at r = r^,. However, the density 
p = F'^ F is ill-defined at r = r^, so it does not make 
sense to consider the continuity condition for jr at this 
point. 

Since the wave functions F± {E, r) are not orthogonal 
to each other, it is convenient to introduce the dual ba- 
sis functions G± {E, r) being the linear combinations of 
F± {E, r) and satisfying the relations: 



00 

/ 



dr rGl {E, r) F,, {E' , r) = 2tt6 {E - E') 5..^ . (58) 



From Eqs. ([531), (EZ]) we have 



G±{E,r) 



F± - iF^ 
1-72 



(59) 



Using Eq. ((5^ . we find that the large-distance asymp- 
totics of G± corresponds to the convergent/divergent 
spherical wave, respectively: 



gT»(B'-+¥')/2 / 



(60) 



This remarkable property is important in the considera- 
tion of the time-dependent problem. 



IV. TIME-DEPENDENT PROBLEM 

The stationary solutions for Af 7^ 0, derived in the 
previous section, look very unusual due to a singular be- 
haviour at r = r^. This behaviour should be reflected in 



the time evolution of wave packets. Note that Eq. p6)) 
has the form resolved with respect to the derivative dtF. 
Besides, the coefficients of the differential operator Hr in 
the right-hand side are smooth functions on the interval 
(0,cx)). So, for suitable initial and boundary conditions, 
the problem of finding the solution F (i, r) is well-posed. 

We choose the initial conditions as 




n{x) = {I ~ x'^f e {I - x'') , 



(62) 



where C is the normalization constant, determined by the 
relation J„ drr |_F| =1. This form corresponds to the 
scattering problem setup and describes the wave packet 
with the average energy E^ and the width A, moving 
from the large distance ro towards the origin (cf. Eq. 
([5^). We assume that ro S> A 3> |-Eo| , i-c, the packet 
width in the coordinate space is small compared to the 
average value of r , and the width in the momentum 
space is small compared to the average value oi E. If the 
potential V (r) is a localized function falling off at r ~ R, 
then we also assume that A ^ i?. Performing the numer- 
ical integration of Eq. (fTH]) over t, we find F{t,r). The 
initial packet with the energy Eq well above Vniax (when 
Eq - Kiax > 1/A) or well below zero (-£'0 > 1/A) 
moves with the speed 2vf (2 in our units), comes to 
small distances, reflects, and goes away. The norm of the 
outgoing packet is the same as that of the incoming one. 
This behaviour looks very similar to that of the wave 
packet obeying the massless Dirac equation in the cen- 
tral external field. The evolution of the wave packet with 
the energy Eq deep inside the interval (0, Knax) is essen- 
tially different. During the scattering process a narrow 
peak develops at r^, (-Bo)- The form of the peak stabilizes 
at large time. The norm of the outgoing packet (corre- 
sponding to the reflected particles) is less than that of 
the incoming one. However, the total norm is conserved 
due to the additional contribution of the peak at finite 
distances (corresponding to the "adhered" particles). To 
demonstrate this behaviour, the time evolution of the 
wave packet (|6T|) in the Coulomb potential V (r) = a/r 
is shown in Fig. [T] One can see that, after refiection, a 
narrow peak appears at r^, = a/ Eq. 

Decomposition method. In order to gain deeper in- 
sight into this behaviour, let us derive the time evolution 
of the wave packet using the decomposition of the ini- 
tial wave packet over the stationary wave functions. The 
decomposition has the form 



r\F[ 
0.5 r 




r\F[ 
0.5 r 




Figure 1: Time evolution of the density, corresponding to the wave packet (|6ip . in the Coulomb potential V (r) = a/r. The 
parameters are , Eo = 1, ro = 60, A = 20,a = 2, M = l (left) and Af = 2 (right). Insets: the form of the peak at large time. 



Fit,r) = 



Vmax 

J £'"'''* [^+ (^^ ^+ (^' '') + ^- (^) ^- (^' ")] 



oo 



(63) 



The coefficients C (E) and C± [E) have the form 



C{E) ^ / dr rFl {E, r) F (0, r) , (64) 



OO 

C±(£;) = I drrGi{E,r)F{0,r). (65) 



Note that the coefficients C± (-E) in front of F± {E, r) 
in the decomposition (j63p are expressed via the overlap 
integrals of F(0,r) vifith the elements G\. {E,r) of the 



r 



dual basis. Let us consider the decomposition of the wave 
packet (|6T|) vifith the energy Eq deep inside the interval 
(0, Vmax), when Kiax-^o > 1/A and Eo > 1/A. In this 
case the main contribution to the integrals in Eqs. (j64p 
and (j65p comes from large distances r '^ tq. Therefore, 
for the calculation of the coefficients C (E) and C± (E), 
we can use the large-r asymptotics (j45|) and (f52|) . We 
obtain that C (ii^) and C_ (£^) are suppressed due to the 
fast oscillations of the integrands, and we can omit the 
corresponding contributions in Eq. (|63p. The coefficient 
C+ (E) has the form 



C+ [E) = C [E) exp {iEro/2) Cl 



{E-Eo)A 



/16 
dxexp (iqx) n{x) = —^ [(3 — g^) sing — 3(7 cos g] 

CiE)-- 



VA72 



-i{Eoro-'p{E))/2 



J dxn^ (a 



315A 
512 ' 



i(£;oi-o-v(-E))/2 



(66) 
(67) 

(68) 



The function f] / (^ 2°^'^ ) ^^ peaked around E = Eq slowly varying function of the energy. Therefore, we can 
with the characteristic width 1/A, while C [E) is some 
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represent the function F {t, r) as 



asymptotic form of F {t, r) for t satisfying the condi- 
tion 1 2t — rg I ^ A and for r obeying the condition 
r* \V' (r*)| \r - r^\ < 1. Keeping in F+ {E,r), Eq. gZl) , 
only the singular component 



F {t, r)=C iEo)J ^e-'-f2 (^^^|^) ^+ (^' ^ ' 



(69) 
where t = t — ro/2. 

Let us demonstrate now that this decomposition leads 
to the appearance of the peak in the vicinity of r = 
Ti, (Eq) at large t. For this purpose we consider the and using the identity (j50p . we obtain 

I 



f+iE,r) 



2Md+ {E, r) 
r{V{r)-E-iO) ' 



(70) 





where V = V [r). Passing to the variable e = E — V we have 

L2 



(71) 



/(,.)« ^^^^^^e- 



de 

2^' 



■iVr I rXg-^^^fl 



Ll 



^"^-/° + -^' )(.,M.)-pf|^.(K + .,.), 



(72) 



where Li = — F, L2 = Knax — 1^- Since |t| ^ A and |tLi^2| ^ 1, we can write Eq. (|72|) as 



/„,,, . H^££M,-.wfi ('A(i_^x ,^ „,„, I *^..„ A ,,^, _ pi 



(73) 



The remaining integral is equal to iO (r), and finally we come to the asymptotics of / (f, r) at \t\ — \t — ?'o/2| ^ A: 

(74) 



/ (, ., . 3!M6(iW ,-.,.0 f ^S;^ I .. (K .) -> (r) 



r 



Note that the right-hand side of Eq. ([74| vanishes for 
T <C —A , and the leading asymptotics of F (i, r) comes 
from the contribution of the nonsingular terms. How- 
ever, we can claim that this asymptotics is not peaked 
in the vicinity of r* {Eq). For r ^ A, the density 
|F| « I/I is independent of r and peaked, due to 
the factor J7 in Eq. (TMl) . with the characteristic width 
5^1/ \V' [Vi,) A|. Thus, we have demonstrated that the 
decomposition method leads to the appearance of the 
peak at large time, which is in agreement with the result 
of direct numerical solution of the differential equation. 
For the Coulomb potential, we have also checked numer- 
ically that the time evolution, obtained by the decompo- 
sition method, coincides with that obtained by the direct 
numerical solution of the differential equation, see Fig. [1] 

Adhesion coefficient. Let us consider the quantity 



lim 

t— >--foo 



Li 



drr\F{t,', 



(75) 



The upper limit L in this formula is some fixed parameter 
obeying the condition L ^ r^, (Eq). The quantity A is 
the adhesion coefficient, i.e., the probability for the two 
particles to remain at finite distances at large t. Using 
Eqs. dZl), dSZl), dSgi, dMl), and dSH) we obtain 



Ak hm 

t— >-|-oo 



drr\fit,r)\' 



(76) 



where 7 = '-f{Eo). For the Coulomb potential V (r) — 
a/r, the dimensional arguments lead to the independence 
of the quantity 7 of |-Eo|- In Fig. [2] the adhesion coeffi- 
cient A for the case of Coulomb potential is shown as a 
function of a. One can see that A grows when a changes 
from to its critical value aM — Vl + 4M2 (when the 
parameter v in Eq. p2p vanishes). 

Asymptotically localized state. The numerical solu- 
tion of the differential equation P^ shows that at large t 
and fixed r the functions h (t, r) and d {t, r) vanish. Then, 

it follows from Eq. ^ that f{t,r) *-=^ e-*^('^)*/o (r). 
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Figure 2: The adhesion coefficient A, Eq. (|76|l. in the 
Coulomb potential as a function of a for M = 1 (solid curve), 
M = 2 (dashed curve), M = 3 (dash-dotted curve). The ver- 
tical lines correspond to the critical values Om = Vl + 4M2. 



where /o (r) is some function of r. This form of the 
asymptotics is also in agreement with Eq. (TMl) . The 
function /q (r) is peaked at r = r^, and depends on the 
form of the initial packet (pT|) . This asymptotics could 
be considered as a hint for the existence of the normal- 
izable solutions of Eq. (IT51) with the constant density 
p(r) = |/o(r)| . However, the equation (IT5|) has no solu- 
tions of the form 



Fit,r) 



^ g-»V(r)t 




Indeed, this form satisfies Eqs. (IT51) and ([H]), but not 
P^ . Instead, we search the asymptotics of the solution 
of Eqs. (HSJ-dUD as 



(77) 



Substituting this form in Eq. (|16l) . we obtain the recur- 
rence relations for /„ (r), ft,„ (r), and d„ (r) which can 
be used to express the asymptotics (|77p via one function 
/o (r). In the leading order we obtain 




f{t,r)^ 


e 


iV{r)t 


h{t,r)^ 


e' 


-iV(r)t 




t 


d{t,r)^ 


e' 


-iV(r)t 



[/o(r) + 0(t-^)] 
iMfo (r) 



rV (r) 
2rt/'2 (r) 






(78) 
(79) 

(80) 



Though the function h {t, r) vanishes at t — > oo , 
its derivative drh(t,r) does not vanish, \drh{t,r)\ ^■ 
\Mf (r) /r\. Due to this behaviour of h {t, r), the equa- 
tion (ITSI) is now satisfied. We see that the asymptotics 



(l78l) - ((80)) is in agreement with the behaviour observed in 
the numerical solution of the differential equation. Note 
that the asymptotics (f78| -(f80 |) leads to the vanishing ra- 
dial and azimuthal current, Eq. ([^0]). 



V. CONCLUSION 

In the present paper, we have considered a model of 
the electron-electron interaction in graphene, based on 
the hamiltonian ([3]). Despite the simplicity of the model, 
it leads to a very unusual dynamics. We have shown, 
both numerically and analytically, that in the process of 
the wave packet scattering the asymptotically localized 
state appears, see Fig. ((TJ and Eq. (TMl) . From the point 
of view of the outside observer, the scattering seems to 
be inelastic. The origin of this phenomenon is traced 
back to two-fold degeneracy of the spectrum of hamilto- 
nian Hr, Eq. ((ni), at < £; < Kiax, Eqs. ^ and 
(l43t . This degeneracy is related to the degeneracy of the 
derivative matrix in Hr- The results obtained are valid 
for any smooth monotonically decreasing vanishing po- 
tential. Though we did not take into account the Fermi 
statistics of the interacting electrons, the requirement of 
the Fermi statistics can be satisfied by appropriate choice 
of the spin part of the wave function. 

For simplicity, we considered the states with definite 
value of J^. It is obvious, that the observed phenomenon 
(the appearance of the asymptotically localized state) 
also retains for any superposition of the states with dif- 
ferent values of J^. Our consideration is not directly 
applicable to the case P 7^ 0, but the appearance of the 
asymptotically localized state is likely to take place also 
in this case because the derivative matrix in Hy, Eq. ([3]), 
is also degenerate. 

We emphasize that this simple model does not take 
into account the existense of the electrons below Fermi 
surface. The effect of such electrons is the appearance 
of virtual electron-hole pairs in the intermediate states. 
While the exact account of this effect is hardly possible, 
its qualitative consideration is very important and will 
be presented elsewhere. It seems that the effect of the 
electrons below Fermi surface is small at least in the case 
when SEp < and \SEf\ » Eq, where SEp is the dif- 
ference between the Fermi energy and the energy of the 
Dirac point. Therefore, it should be possible to observe 
the appearance of the asymptotically localized state in 
the experiment. Note that the differential equation ([5]) 
does not take into account the effect of the finite lattice 
scale lo ^ 0.14nni, and valid if the wave function changes 
slowly on this scale. It gives us two conditions: the width 
of the localized state should be much larger than the lat- 
tice scale, and the variation of the phase on the lattice 
scale should be small compared to unity. The first con- 
dition gives the constraint on the energy dispersion SE 
in the wave packet: 



SE:^\V'{r,)\lo. 
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It follows from Eq. ((74)) that, at large time, the phase of 
the wave function varies significantly on the lattice scale. 
The second condition gives us the lower bound tq of the 
life time of the localized state: 



1 



To 



^o|^'(r.)| 



>l/<5^. 



which means that the localized state lives long enough. 
For the Coulomb potential, these two conditions read 



SE > 



To 



a 
a 



Note that the first condition is compatible with the con- 
dition SE <^ Eq provided that Eq is sufficiently small. 
For instance, for Eq ~ ImeV and a ~ 1, we have quite 
large time tq ~ 1/is. Therefore, if the model considered is 
relevant to the electron-electron interaction in graphene, 
one may hope to observe the long-lived localized states 
in the experiment. 
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